#include <fstream>
#include <iostream>
#include <sstream>
#include <string>

#include "utils.h"
#include "GRDECLData.h"

// Print the values of the Corn8 object
void Corn8::print() {
    std::cout << zval[0][0][0] << "--" << zval[0][0][1] << "\t\t" << zval[1][0][0] << "--" << zval[1][0][1] << std::endl;
    std::cout << zval[0][1][0] << "--" << zval[0][1][1] << "\t\t" << zval[1][1][0] << "--" << zval[1][1][1] << std::endl;
}

// Calculate and return the center point of the Corn8 object
Point3d Corn8::GetCenter() {
    Point3d res(0, 0, 0);
    for(int l = 0; l < 2; l++) {
        for(int m = 0; m < 2; m++) {
            for(int n = 0; n < 2; n++) {
                res += coord[l][m][n];
            }
        }
    }
    return res / 8;
}

// Print the properties of the grid
void GridProperties::print() {
    std::cout << "Grid Properties:" << std::endl;
    std::cout << "Pillar Vertical: " << isPillarVertical << std::endl;
    std::cout << "Grid Dense: " << isGridDense << std::endl;
    std::cout << "Coord Monotonic: " << isCoordMonotonic << std::endl;
}

// Print the range of the grid
void GridRange::print() {
    std::cout << "Grid Range:" << std::endl;
    std::cout << "rmin: " << rmin.x << ", " << rmin.y << ", " << rmin.z << std::endl;
    std::cout << "rmax: " << rmax.x << ", " << rmax.y << ", " << rmax.z << std::endl;
}

Line& GRDECLData::GetPillar(INDEX y, INDEX x) {
    return pillars[y * nx1 + x];
}

// Get the corner at the specified (z, y, x) index
Corn8& GRDECLData::GetCorner(INDEX z, INDEX y, INDEX x) {
    return corns[z * nxy + y * nx + x];
}

// Get the Young's modulus at the specified (z, y, x) index
double GRDECLData::GetYoung(INDEX z, INDEX y, INDEX x) {
    return young[z * nxy + y * nx + x];
}

// Get the Poisson's ratio at the specified (z, y, x) index
double GRDECLData::GetPoisson(INDEX z, INDEX y, INDEX x) {
    return poisson[z * nxy + y * nx + x];
}

double GRDECLData::GetRockDen(INDEX z, INDEX y, INDEX x) {
    return rockden[z * nxy + y * nx + x];
}

double GRDECLData::GetPoro(INDEX z, INDEX y, INDEX x) {
    return poro[z * nxy + y * nx + x];
}

// Constructor that initializes grid data from a file
GRDECLData::GRDECLData(std::string filename) : filename(filename) {
    std::cout << "Init grid data from file " << filename << std::endl;
    if(InitFromFileOpt()) {
        std::cerr << "Failed to init from file: " << filename << std::endl;
        // quit program
        exit(-1);
    }
    InitProperties();
    InitRange();
    // InitConvex();
}

geo2d::Point2d GRDECLData::GetPillarAtZval(INDEX y, INDEX x, double zval) {
    geo2d::Point2d res;
    Line line = GetPillar(y, x);
    double t = (zval - line.a.z) / (line.b.z - line.a.z);
    res.x = line.a.x + t * (line.b.x - line.a.x);
    res.y = line.a.y + t * (line.b.y - line.a.y);
    return res;
}

geo2d::Point2d GRDECLData::GetPillarAtZval(const Line& pillar, double zval) {
    geo2d::Point2d res;
    double t = (zval - pillar.a.z) / (pillar.b.z - pillar.a.z);
    res.x = pillar.a.x + t * (pillar.b.x - pillar.a.x);
    res.y = pillar.a.y + t * (pillar.b.y - pillar.a.y);
    return res;
}

Point2d GRDECLData::GetPointFromPillar(INDEX y, INDEX x, double zval) {
    Point2d res;
    Line line = GetPillar(y, x);
    double t = (zval - line.a.z) / (line.b.z - line.a.z);
    res.x = line.a.x + t * (line.b.x - line.a.x);
    res.y = line.a.y + t * (line.b.y - line.a.y);
    return res;
}

/**
* Optimized version of InitFromFile using IOBuffer.
* return 0 if success, -1 if failed
*/
int GRDECLData::InitFromFileOpt() {
    IOBuffer buffer(filename);

    // Check if the file is finished or cannot be opened
    if(buffer.file_finished) {
        std::cerr << "Unable to open file: " << filename << std::endl;
        return -1;
    }

    std::string sec_name; // section name

    // Loop through each section in the file
    while(buffer.getString(sec_name)) {
        // Check for MAPUNITS section
        if(sec_name.find("MAPUNITS") != std::string::npos) {
            std::string unit_name;
            buffer.getString(unit_name);
            // Ensure the units are in meters
            if(unit_name.find("METRES") == std::string::npos) {
                std::cerr << unit_name << std::endl;
                std::cerr << "Map units not meters" << std::endl;
                return -1;
            }
        }

        if (sec_name.find("MAPAXES") != std::string::npos) {
            double x1, y1, x2, y2, x3, y3;
            buffer.getFloat(x1);
            buffer.getFloat(y1);
            buffer.getFloat(x2);
            buffer.getFloat(y2);
            buffer.getFloat(x3);
            buffer.getFloat(y3);
            origin_x = x2;
            origin_y = y2;
            printf("origin_x = %f, origin_y = %f\n", origin_x, origin_y);
        }

        // Check for SPECGRID section
        if(sec_name.find("SPECGRID") != std::string::npos) {
            buffer.getInteger(nx);
            buffer.getInteger(ny);
            buffer.getInteger(nz);
            printf("nx=%d, ny=%d, nz=%d\n", nx, ny, nz);
            nxy = nx * ny;
            nxyz = nx * ny * nz;
            nx1 = nx + 1;
            // Allocate memory for grid data
            pillars = new Line[(nx + 1) * (ny + 1)];
            corns = new Corn8[nxyz];
            young = new double[nxyz];
            poisson = new double[nxyz];
            rockden = new double[nxyz];
            poro = new double[nxyz];
            actNum = new bool[nxyz];
            isDense = new bool[nxyz];
            isFault = new bool[(nx + 1) * (ny + 1) * (nz + 1)];
            isConvex = new bool[nx * ny];

            for (int i = 0; i < nxyz; i++) {
                isDense[i] = true;
            }
            for (int i = 0; i < (nx + 1) * (ny + 1) * (nz + 1); i++) {
                isFault[i] = false;
            }
            for (int i = 0; i < nx * ny; i++) {
                isConvex[i] = true;
            }
        }

        // Check for COORD section
        if(sec_name.find("COORD") != std::string::npos) {
            // Read coordinates for each pillar
            for(INDEX j = 0; j < ny + 1; j++) {
                for(INDEX i = 0; i < nx + 1; i++) {
                    buffer.getFloat(GetPillar(j, i).a.x);
                    buffer.getFloat(GetPillar(j, i).a.y);
                    buffer.getFloat(GetPillar(j, i).a.z);
                    // printf("%f %f %f\n", GetPillar(j, i).a.x, GetPillar(j, i).a.y, GetPillar(j, i).a.z);

                    buffer.getFloat(GetPillar(j, i).b.x);
                    buffer.getFloat(GetPillar(j, i).b.y);
                    buffer.getFloat(GetPillar(j, i).b.z);
                    // printf("%f %f %f\n", GetPillar(j, i).b.x, GetPillar(j, i).b.y, GetPillar(j, i).b.z);
                }
            }
        }

        // Check for ZCORN section
        if(sec_name.find("ZCORN") != std::string::npos) {
            int cnt = 0;
            // Read Z coordinates for each corner
            for(INDEX k = 0; k < nz * 2; k++) {
                for(INDEX j = 0; j < ny * 2; j++) {
                    for(INDEX i = 0; i < nx * 2; i++) {
                        double val;
                        buffer.getFloat(val);
                        // printf("%.10lf ", val);
                        // std::cout << val << "\t";
                        // if(++cnt % 6 == 0) {
                        //     printf("\n");
                        // }
                        // fflush(stdout);
                        GetCorner(k/2, j/2, i/2).zval[k & 1][j & 1][i & 1] = val;
                    }
                }
            }
            // Calculate coordinates for each corner
            for(INDEX k = 0; k < nz; k++) {
                for(INDEX j = 0; j < ny; j++) {
                    for(INDEX i = 0; i < nx; i++) {
                        for(int l = 0; l < 2; l++) {
                            for(int m = 0; m < 2; m++) {
                                for(int n = 0; n < 2; n++) {
                                    double portion = (GetCorner(k, j, i).zval[l][m][n] - GetPillar(j + m, i + n).a.z) / (GetPillar(j + m, i + n).b.z - GetPillar(j + m, i + n).a.z);
                                    GetCorner(k, j, i).coord[l][m][n].x = GetPillar(j + m, i + n).a.x + portion * (GetPillar(j + m, i + n).b.x - GetPillar(j + m, i + n).a.x);
                                    GetCorner(k, j, i).coord[l][m][n].y = GetPillar(j + m, i + n).a.y + portion * (GetPillar(j + m, i + n).b.y - GetPillar(j + m, i + n).a.y);
                                    GetCorner(k, j, i).coord[l][m][n].z = GetCorner(k, j, i).zval[l][m][n];
                                }
                            }
                        }
                    }
                }
            }
        }

        std::string line;
        auto ReadProp = [&](bool* prop) {
            INDEX cnt = 0;
            while(cnt < nxyz) {
                int num_val;
                double val;
                if(!buffer.getExp(num_val, val)) {
                    break;
                }
                while(num_val) {
                    prop[cnt++] = val;
                    --num_val;
                }
            }
        };

        if (sec_name.find("ACTNUM") != std::string::npos) {
            ReadProp(actNum);
        }
    }
    return 0;
}

// Initialize the properties of the GRDECLData object void
void GRDECLData::InitProperties() {
    // Set initial properties 
    properties.inited = true;
    properties.isPillarVertical = true;
    properties.isGridDense = true;
    properties.isCoordMonotonic = true;

    // Check if all pillars are vertical
    for(int j = 0; j < ny + 1; j++) {
        for(int i = 0; i < nx + 1; i++) {
            if(GetPillar(j, i).a.x != GetPillar(j, i).b.x || GetPillar(j, i).a.y != GetPillar(j, i).b.y) {
                properties.isPillarVertical = false;
            }
        }
    }
    
    int dense_cnt = 0;
    // Check if the grid is dense
    for(INDEX k = 0; k < nz; k++) {
        for(INDEX j = 0; j < ny; j++) {
            for(INDEX i = 0; i < nx; i++) {
                // Check density in x direction
                if(i + 1 < nx) {
                    if (!isEqual(GetCorner(k, j, i).zval[0][0][1], GetCorner(k, j, i+1).zval[0][0][0]) ||
                        !isEqual(GetCorner(k, j, i).zval[1][0][1], GetCorner(k, j, i+1).zval[1][0][0]) ||
                        !isEqual(GetCorner(k, j, i).zval[0][1][1], GetCorner(k, j, i+1).zval[0][1][0]) ||
                        !isEqual(GetCorner(k, j, i).zval[1][1][1], GetCorner(k, j, i+1).zval[1][1][0])) {
                        if(properties.isGridDense) {
                            // printf("x not dense at (%d, %d, %d) with (%d, %d, %d)\n", k, j, i, k, j, i+1);
                            // GetCorner(k, j, i).print();
                            // GetCorner(k, j, i+1).print();
                            dense_cnt++;
                            isDense[k * nxy + j * nx + i] = false;
                            isDense[k * nxy + j * nx + i + 1] = false;
                            // continue;
                            // properties.isGridDense = false;
                        }
                    }
                }
                // Check density in y direction
                if(j + 1 < ny) {
                    if (!isEqual(GetCorner(k, j, i).zval[0][1][0], GetCorner(k, j+1, i).zval[0][0][0]) ||
                        !isEqual(GetCorner(k, j, i).zval[1][1][0], GetCorner(k, j+1, i).zval[1][0][0]) ||
                        !isEqual(GetCorner(k, j, i).zval[0][1][1], GetCorner(k, j+1, i).zval[0][0][1]) ||
                        !isEqual(GetCorner(k, j, i).zval[1][1][1], GetCorner(k, j+1, i).zval[1][0][1])) {
                        if(properties.isGridDense) {
                            // printf("y not dense at (%d, %d, %d) with (%d, %d, %d)\n", k, j, i, k, j+1, i);
                            // GetCorner(k, j, i).print();
                            // GetCorner(k, j+1, i).print();
                            dense_cnt++;
                            isDense[k * nxy + j * nx + i] = false;
                            isDense[k * nxy + (j + 1) * nx + i] = false;
                            // continue;
                            // properties.isGridDense = false;
                        }
                    }
                }
                // Check density in z direction
                if(k + 1 < nz) {
                    if (!isEqual(GetCorner(k, j, i).zval[1][0][0], GetCorner(k+1, j, i).zval[0][0][0]) ||
                        !isEqual(GetCorner(k, j, i).zval[1][1][0], GetCorner(k+1, j, i).zval[0][1][0]) ||
                        !isEqual(GetCorner(k, j, i).zval[1][0][1], GetCorner(k+1, j, i).zval[0][0][1]) ||
                        !isEqual(GetCorner(k, j, i).zval[1][1][1], GetCorner(k+1, j, i).zval[0][1][1])) {
                        if(properties.isGridDense) {
                            // printf("z not dense at (%d, %d, %d) with (%d, %d, %d)\n", k, j, i, k+1, j, i);
                            // GetCorner(k, j, i).print();
                            // GetCorner(k+1, j, i).print();
                            dense_cnt++;
                            isDense[k * nxy + j * nx + i] = false;
                            isDense[(k + 1) * nxy + j * nx + i] = false;
                            // continue;
                            // properties.isGridDense = false;
                        }
                    }
                }
            }
        }
    }

    printf("dense count: %d\n", dense_cnt);

    // Determine monotonic direction for x, y, and z coordinates
    properties.x_mono_dir = GetCorner(0, 0, 0).coord[0][0][0].x < GetCorner(0, 0, 0).coord[0][0][1].x ? 1 : -1;
    properties.y_mono_dir = GetCorner(0, 0, 0).coord[0][0][0].y < GetCorner(0, 0, 0).coord[0][1][0].y ? 1 : -1;
    properties.z_mono_dir = GetCorner(0, 0, 0).coord[0][0][0].z < GetCorner(0, 0, 0).coord[1][0][0].z ? 1 : -1;

    // Check if coordinates are monotonic
    for(INDEX k = 0; k < nz; k++) {
        for(INDEX j = 0; j < ny; j++) {
            for(INDEX i = 0; i < nx; i++) {
                Corn8 corn = GetCorner(k, j, i);
                if(corn.zval[0][0][0] * properties.z_mono_dir > corn.zval[1][0][0] * properties.z_mono_dir || 
                   corn.zval[0][0][1] * properties.z_mono_dir > corn.zval[1][0][1] * properties.z_mono_dir || 
                   corn.zval[0][1][0] * properties.z_mono_dir > corn.zval[1][1][0] * properties.z_mono_dir || 
                   corn.zval[0][1][1] * properties.z_mono_dir > corn.zval[1][1][1] * properties.z_mono_dir) {
                    // std::cout << "zval not monotonic at (" << k << ", " << j << ", " << i << ")" << std::endl;
                    // std::cout << corn.zval[0][0][0] << "--" << corn.zval[1][0][0] << "\t\t" << corn.zval[0][0][1] << "--" << corn.zval[1][0][1] << std::endl;
                    // std::cout << corn.zval[0][1][0] << "--" << corn.zval[1][1][0] << "\t\t" << corn.zval[0][1][1] << "--" << corn.zval[1][1][1] << std::endl;
                    properties.isCoordMonotonic = false;
                }
            }
        }
    }
    for(int j = 0; j < ny; j++) {
        for(int i = 0; i < nx; i++) {
            if(GetPillar(j, i).a.x * properties.x_mono_dir > GetPillar(j, i + 1).a.x * properties.x_mono_dir || 
               GetPillar(j, i).a.y * properties.y_mono_dir > GetPillar(j + 1, i).a.y * properties.y_mono_dir ||
               GetPillar(j, i).b.x * properties.x_mono_dir > GetPillar(j, i + 1).b.x * properties.x_mono_dir || 
               GetPillar(j, i).b.y * properties.y_mono_dir > GetPillar(j + 1, i).b.y * properties.y_mono_dir ) {
                // std::cout << "pillar not monotonic at (" << j << ", " << i << ")" << std::endl;
                // std::cout << GetPillar(j, i).a.x << "--" << GetPillar(j, i + 1).a.x << "\t\t" << GetPillar(j, i).a.y << "--" << GetPillar(j + 1, i).a.y << std::endl;
                properties.isCoordMonotonic = false;
            }
        }
    }

    properties.print();
}

void GRDECLData::InitRange() {
    // Initialize the minimum range values to a very large number
    range.rmin.x = range.rmin.y = range.rmin.z = 1e19;
    // Initialize the maximum range values to a very small number
    range.rmax.x = range.rmax.y = range.rmax.z = -1e19;

    // In fact, we can only use 2 points
    range.rmin.x = std::min(GetCorner(0, 0, 0).coord[0][0][0].x, GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].x);
    range.rmin.y = std::min(GetCorner(0, 0, 0).coord[0][0][0].y, GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].y);
    range.rmin.z = std::min(GetCorner(0, 0, 0).coord[0][0][0].z, GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].z);
    range.rmax.x = std::max(GetCorner(0, 0, 0).coord[0][0][0].x, GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].x);
    range.rmax.y = std::max(GetCorner(0, 0, 0).coord[0][0][0].y, GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].y);
    range.rmax.z = std::max(GetCorner(0, 0, 0).coord[0][0][0].z, GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].z);

    // Print the calculated range values
    range.print();
}

// Initialize the isConvex array
void GRDECLData::InitConvex() {
    for (INDEX j = 0; j < ny; j++) {
        for (INDEX i = 0; i < nx; i++) {
            Point2d p1 = {GetCorner(0, j, i).coord[0][0][0].x, GetCorner(0, j, i).coord[0][0][0].y};
            Point2d p2 = {GetCorner(0, j, i).coord[0][0][1].x, GetCorner(0, j, i).coord[0][0][1].y};
            Point2d p3 = {GetCorner(0, j, i).coord[0][1][0].x, GetCorner(0, j, i).coord[0][1][0].y};
            Point2d p4 = {GetCorner(0, j, i).coord[0][1][1].x, GetCorner(0, j, i).coord[0][1][1].y};
            auto convex1 = geo2d::checkQuadrilateralConvexity(p1, p2, p3, p4);

            Point2d p5 = {GetCorner(nz - 1, j, i).coord[1][0][0].x, GetCorner(nz - 1, j, i).coord[1][0][0].y};
            Point2d p6 = {GetCorner(nz - 1, j, i).coord[1][0][1].x, GetCorner(nz - 1, j, i).coord[1][0][1].y};
            Point2d p7 = {GetCorner(nz - 1, j, i).coord[1][1][0].x, GetCorner(nz - 1, j, i).coord[1][1][0].y};
            Point2d p8 = {GetCorner(nz - 1, j, i).coord[1][1][1].x, GetCorner(nz - 1, j, i).coord[1][1][1].y};
            auto convex2 = geo2d::checkQuadrilateralConvexity(p5, p6, p7, p8);

            isConvex[j * nx + i] = convex1.first && convex2.first;
            if (!isConvex[j * nx + i]) {
                printf("convexity not satisfied at (%d, %d)\n", j, i);
            }
        }
    }
}

void GRDECLData::ReadPropsFromFileOpt(std::string filename) {
    // Initialize IOBuffer with the given filename
    IOBuffer buffer(filename);
    if(buffer.file_finished) {
        // If the file cannot be opened, print an error message and return
        std::cerr << "Unable to open file: " << filename << std::endl;
        return;
    }

    std::string line;
    // Lambda function to read property values from the buffer
    auto ReadProp = [&](double* prop) {
        INDEX cnt = 0;
        while(cnt < nxyz) {
            int num_val;
            double val;
            // Get the next expression from the buffer
            if(!buffer.getExp(num_val, val)) {
                break;
            }
            // Assign the value to the property array
            while(num_val) {
                prop[cnt++] = val;
                --num_val;
            }
        }
    };

    std::string sec_name; // section name
    // Read section names from the buffer
    while(buffer.getString(sec_name)) {
        // Check for the "YOUNGMOD" keyword and read the Young's modulus values
        if(sec_name.find("YOUNGMOD") != std::string::npos) {
            ReadProp(young);
        }
        // Check for the "POISSONR" keyword and read the Poisson's ratio values
        if(sec_name.find("POISSONR") != std::string::npos) {
            ReadProp(poisson);
        }
        // Check for the "ROCKDEN" keyword and read the rock dense values
        if(sec_name.find("ROCKDEN") != std::string::npos) {
            ReadProp(rockden);
        }
        // Check for the "PORO" keyword and read the poro values
        if(sec_name.find("PORO") != std::string::npos) {
            ReadProp(poro);
        }
    }

    // Calculate the average property values
    CalAverageProps();
}

void GRDECLData::CalAverageProps() {
    // Initialize the average property values
    avg_young = 0;
    avg_poisson = 0;
    avg_rockden = 0;
    avg_poro = 0;

    // Calculate the sum of all property values
    for(INDEX i = 0; i < nxyz; i++) {
        avg_young += young[i];
        avg_poisson += poisson[i];
        avg_rockden += rockden[i];
        avg_poro += poro[i];
    }

    // Calculate the average property values
    avg_young /= nxyz;
    avg_poisson /= nxyz;
    avg_rockden /= nxyz;
    avg_poro /= nxyz;
}

void GRDECLData::PrintHeadProps() {
    // Print the first 10 values of the Young's modulus array
    for(int i = 0; i < 10; i++) {
        printf("young[%d]=%f\n", i, young[i]);
    }
    // Print the first 10 values of the Poisson's ratio array
    for(int i = 0; i < 10; i++) {
        printf("poisson[%d]=%f\n", i, poisson[i]);
    }
    // Print the first 10 values of the rockden array
    for(int i = 0; i < 10; i++) {
        printf("rockden[%d]=%f\n", i, rockden[i]);
    }
    // Print the first 10 values of the poro array
    for(int i = 0; i < 10; i++) {
        printf("poro[%d]=%f\n", i, poro[i]);
    }
}

int GRDECLData::PrintPropToFile(std::string filename) {
    // Open the file for writing
    FILE* file = fopen(filename.c_str(), "w");
    if(file == NULL) {
        // If the file cannot be opened, print an error message and return -1
        std::cerr << "Unable to open file: " << filename << std::endl;
        return -1;
    }

    // Lambda function to write property values to the file
    auto WriteProp = [&](double* prop) {
        INDEX start_pos = 0;
        INDEX end_pos = 0;
        // Loop through all property values
        while(end_pos < nxyz) {
            // Find the range of identical values
            while(end_pos < nxyz && prop[start_pos] == prop[end_pos]) {
                ++end_pos;
            }
            // Write the range of identical values to the file
            if(end_pos - start_pos > 0) {
                fprintf(file, "%d*%f\n", end_pos - start_pos, prop[start_pos]);
            }
            start_pos = end_pos;
        }
    };
    // Write Young's modulus values to the file
    fprintf(file, "YOUNGMOD\n");
    WriteProp(young);
    // Write Poisson's ratio values to the file
    fprintf(file, "POISSONR\n");
    WriteProp(poisson);
    // Close the file
    fclose(file);
    return 0;
}

void GRDECLData::ReadFaultFromFile(std::string filename) {
    // Open the file for reading
    IOBuffer buffer(filename);
    if(buffer.file_finished) {
        // If the file cannot be opened, print an error message and return
        std::cerr << "Unable to open file: " << filename << std::endl;
        return;
    }

    std::string temp_fault, fault_name; // section name
    // Read section names from the buffer
    while(buffer.getString(temp_fault)) {
        buffer.getString(fault_name);
        double t1, t2, t3;
        int v1;

        geo3d::Fault fault;
        fault.name = fault_name;
        // fault.lines.resize(64);
        while(true) {
            char c = buffer.getFloat(t1);
            if(c == 'F' || !c) { // get FAULT
                break;
            }
            buffer.getFloat(t2);
            buffer.getInteger(v1);
            buffer.getFloat(t3);
            // printf("get fault %s line : %f %f %d %f\n", fault_name.c_str(), t1, t2, v1, t3);
            if (v1 > fault.lines.size()) {
                fault.lines.resize(fault.lines.size() + 1);
            }
            fault.lines[v1 - 1].push_back({t1 - origin_x, -(t2 - origin_y), -t3}); // reverse y and z
        }
        fault.lines.shrink_to_fit();
        // printf("fault %s has %d lines\n", fault_name.c_str(), fault.lines.size());
        fault.horizontal_size = fault.lines.size();
        fault.vertical_size = fault.lines[0].size();
        faults.push_back(fault);
    }
}

void GRDECLData::WriteVTK(std::string filename, bool cellData) {
    // Open the VTK file for writing
    std::ofstream vtkFile(filename);
    if (!vtkFile.is_open()) {
        // If the file cannot be opened, print an error message and return
        std::cerr << "Unable to open file: " << filename << std::endl;
        return;
    }

    // Write the VTK file header
    vtkFile << "# vtk DataFile Version 2.0\n";
    vtkFile << "Grid Legacy VTK File\n";
    vtkFile << "ASCII\n";
    vtkFile << "DATASET UNSTRUCTURED_GRID\n";
    vtkFile << "\n";

    // Write the grid points to the VTK file
    vtkFile << "POINTS " << 8 * nx * ny * nz << " float\n";
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                for (int l = 0; l < 2; l++) {
                    for (int m = 0; m < 2; m++) {
                        for (int n = 0; n < 2; n++) {
                            Point3d p = GetCorner(k, j, i).coord[l][m][n];
                            vtkFile << p.x << " " << p.y << " " << p.z << "\n";
                        }
                    }
                }
            }
        }
    }
    vtkFile << "\n";

    // Write the grid cells to the VTK file
    vtkFile << "CELLS " << nx * ny * nz << " " << nx * ny * nz * 9 << "\n";
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << "8 ";
                vtkFile << 8 * (k * ny * nx + j * nx + i) << " "; // cell global index
                vtkFile << 8 * (k * ny * nx + j * nx + i) + 1 << " ";
                vtkFile << 8 * (k * ny * nx + j * nx + i) + 2 << " ";
                vtkFile << 8 * (k * ny * nx + j * nx + i) + 3 << " ";
                vtkFile << 8 * (k * ny * nx + j * nx + i) + 4 << " ";
                vtkFile << 8 * (k * ny * nx + j * nx + i) + 5 << " ";
                vtkFile << 8 * (k * ny * nx + j * nx + i) + 6 << " ";
                vtkFile << 8 * (k * ny * nx + j * nx + i) + 7 << "\n";
            }
        }
    }
    vtkFile << "\n";

    // Write the cell types to the VTK file
    vtkFile << "CELL_TYPES " << nx * ny * nz << "\n"; 
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << "11\n";
            }
        }
    }

    if (!cellData) {
        // Close the VTK file
        vtkFile.close();
        return;
    }

    // Write the cell data to the VTK file
    vtkFile << "CELL_DATA " << nx * ny * nz << "\n";
    vtkFile << "SCALARS young float 1\n";
    vtkFile << "LOOKUP_TABLE default\n";
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << GetYoung(k, j, i) << "\n";
            }
        }
    }
    vtkFile << "\n";

    vtkFile << "SCALARS poisson float 1\n";
    vtkFile << "LOOKUP_TABLE default\n";
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << GetPoisson(k, j, i) << "\n";
            }
        }
    }
    vtkFile << "\n";

    vtkFile << "SCALARS rockden float 1\n";
    vtkFile << "LOOKUP_TABLE default\n";
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << rockden[k * nxy + j * nx + i] << "\n";
            }
        }
    }
    vtkFile << "\n";

    vtkFile << "SCALARS poro float 1\n";
    vtkFile << "LOOKUP_TABLE default\n";
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << poro[k * nxy + j * nx + i] << "\n";
            }
        }
    }
    vtkFile << "\n";

    // vtkFile << "SCALARS actNum int 1\n";
    // vtkFile << "LOOKUP_TABLE default\n";
    // for (INDEX k = 0; k < nz; k++) {
    //     for (INDEX j = 0; j < ny; j++) {
    //         for (INDEX i = 0; i < nx; i++) {
    //             vtkFile << actNum[k * nxy + j * nx + i] << "\n";
    //         }
    //     }
    // }
    // vtkFile << "\n";

    // vtkFile << "SCALARS isDense int 1\n";
    // vtkFile << "LOOKUP_TABLE default\n";
    // for (INDEX k = 0; k < nz; k++) {
    //     for (INDEX j = 0; j < ny; j++) {
    //         for (INDEX i = 0; i < nx; i++) {
    //             vtkFile << isDense[k * nxy + j * nx + i] << "\n";
    //         }
    //     }
    // }
    // vtkFile << "\n";

    // Close the VTK file
    vtkFile.close();
}

void GRDECLData::WriteLayerVTK(std::string base_filename, bool cellData) {
    
    for(INDEX k = 0; k < nz; k++) {
        std::string filename = base_filename + std::to_string(k) + ".vtk";
        // Open the VTK file for writing
        std::ofstream vtkFile(filename);
        if (!vtkFile.is_open()) {
            // If the file cannot be opened, print an error message and continue
            std::cerr << "Unable to open file: " << filename << std::endl;
            continue;
        }

        // Write the VTK file header
        vtkFile << "# vtk DataFile Version 2.0\n";
        vtkFile << "Grid Legacy VTK File\n";
        vtkFile << "ASCII\n";
        vtkFile << "DATASET UNSTRUCTURED_GRID\n";
        vtkFile << "\n";

        // Write the grid points to the VTK file
        vtkFile << "POINTS " << 8 * nx * ny << " float\n";
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                for (int l = 0; l < 2; l++) {
                    for (int m = 0; m < 2; m++) {
                        for (int n = 0; n < 2; n++) {
                            Point3d p = GetCorner(k, j, i).coord[l][m][n];
                            vtkFile << p.x << " " << p.y << " " << p.z << "\n";
                        }
                    }
                }
            }
        }
        vtkFile << "\n";

        // Write the grid cells to the VTK file
        vtkFile << "CELLS " << nx * ny << " " << nx * ny * 9 << "\n";
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << "8 ";
                vtkFile << 8 * (j * nx + i) << " "; // cell global index
                vtkFile << 8 * (j * nx + i) + 1 << " ";
                vtkFile << 8 * (j * nx + i) + 2 << " ";
                vtkFile << 8 * (j * nx + i) + 3 << " ";
                vtkFile << 8 * (j * nx + i) + 4 << " ";
                vtkFile << 8 * (j * nx + i) + 5 << " ";
                vtkFile << 8 * (j * nx + i) + 6 << " ";
                vtkFile << 8 * (j * nx + i) + 7 << "\n";
            }
        }
        vtkFile << "\n";

        // Write the cell types to the VTK file
        vtkFile << "CELL_TYPES " << nx * ny << "\n"; 
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << "11\n";
            }
        }

        if (!cellData) {
            // Close the VTK file
            vtkFile.close();
            continue;
        }

        // Write the cell data to the VTK file
        vtkFile << "CELL_DATA " << nx * ny << "\n";
        vtkFile << "SCALARS young float 1\n";
        vtkFile << "LOOKUP_TABLE default\n";
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << GetYoung(k, j, i) << "\n";
            }
        }
        vtkFile << "\n";

        vtkFile << "SCALARS poisson float 1\n";
        vtkFile << "LOOKUP_TABLE default\n";
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << GetPoisson(k, j, i) << "\n";
            }
        }
        vtkFile << "\n";

        vtkFile << "SCALARS rockden float 1\n";
        vtkFile << "LOOKUP_TABLE default\n";
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << rockden[k * nxy + j * nx + i] << "\n";
            }
        }
        vtkFile << "\n";

        vtkFile << "SCALARS poro float 1\n";
        vtkFile << "LOOKUP_TABLE default\n";
        for (INDEX j = 0; j < ny; j++) {
            for (INDEX i = 0; i < nx; i++) {
                vtkFile << poro[k * nxy + j * nx + i] << "\n";
            }
        }
        vtkFile << "\n";

        // Close the VTK file
        vtkFile.close();        
    }
}

void GRDECLData::WriteFaultVTK(std::string filename) {
    // Open the VTK file for writing
    std::ofstream vtkFile(filename);
    if (!vtkFile.is_open()) {
        // If the file cannot be opened, print an error message and return
        std::cerr << "Unable to open file: " << filename << std::endl;
        return;
    }

    // Write the VTK file header
    vtkFile << "# vtk DataFile Version 2.0\n";
    vtkFile << "Grid Legacy VTK File\n";
    vtkFile << "ASCII\n";
    vtkFile << "DATASET UNSTRUCTURED_GRID\n";
    vtkFile << "\n";

    int points_count = 0;
    for (auto& fault : faults) {
        for (auto& line : fault.lines) {
            points_count += line.size();
        }
    }

    // Write the grid points to the VTK file
    vtkFile << "POINTS " << points_count << " float\n";
    for (auto& fault : faults) {
        for (auto& line : fault.lines) {
            for (auto& point : line) {
                vtkFile << point.x << " " << point.y << " " << point.z << "\n"; // reverse y and z before
            }
        }
    }
    vtkFile << "\n";

    int surfs_count = 0;
    for (auto& fault : faults) {
        int inner_surfs = fault.lines[0].size() - 1;
        surfs_count += (fault.lines.size() - 1) * inner_surfs;
    }

    int cells_count = surfs_count;

    // surface data
    vtkFile << "CELLS " << cells_count << " " << cells_count * 5 << "\n"; 

    int cell_offset = 0;
    for (auto& fault : faults) {
        for (int i = 1; i < fault.lines.size(); i++) {
            for (int j = 0; j < fault.lines[i].size() - 1; j++) {
                vtkFile << "4 ";
                vtkFile << cell_offset + j << " ";
                vtkFile << cell_offset + j + 1 << " ";
                vtkFile << cell_offset + j + fault.lines[i].size() + 1 << " ";
                vtkFile << cell_offset + j + fault.lines[i].size() << "\n";
            }
            cell_offset += fault.lines[i].size();
        }
        cell_offset += fault.lines[0].size();
    }
    vtkFile << "\n";

    // Write the cell types to the VTK file
    vtkFile << "CELL_TYPES " << cells_count << "\n";
    for (int i = 0; i < cells_count; i++) {
        vtkFile << "9\n";
    }
    vtkFile << "\n";

    // Close the VTK file
    vtkFile.close();
}

void GRDECLData::InterpolateNewPillar(double x, double y, double z, Line& new_pillar, INDEX& px, INDEX& py) {
    // 找到 x, y, z 相邻的 pillar
    geo2d::Point2d temp_pillar, temp_pillar1;
    if(_curr_pillar_x + 1 < nx)
        temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x + 1, z);
    while(_curr_pillar_x + 1 < nx && temp_pillar.x * properties.x_mono_dir <= x * properties.x_mono_dir) {
        _curr_pillar_x++;
        temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x + 1, z);
    }
    temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x, z);
    while(_curr_pillar_x - 1 >= 0 && temp_pillar.x * properties.x_mono_dir > x * properties.x_mono_dir) {
        _curr_pillar_x--;
        temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x, z);
    }
    if(_curr_pillar_y + 1 < ny)
        temp_pillar = GetPillarAtZval(_curr_pillar_y + 1, _curr_pillar_x, z);
    while(_curr_pillar_y + 1 < ny && temp_pillar.y * properties.y_mono_dir <= y * properties.y_mono_dir) {
        _curr_pillar_y++;
        temp_pillar = GetPillarAtZval(_curr_pillar_y + 1, _curr_pillar_x, z);
    }
    temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x, z);
    while(_curr_pillar_y - 1 >= 0 && temp_pillar.y * properties.y_mono_dir > y * properties.y_mono_dir) {
        _curr_pillar_y--;
        temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x, z);
    }
    // 重复一遍，防止单维坐标局部问题
    if(_curr_pillar_x + 1 < nx)
        temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x + 1, z);
    while(_curr_pillar_x + 1 < nx && temp_pillar.x * properties.x_mono_dir <= x * properties.x_mono_dir) {
        _curr_pillar_x++;
        temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x + 1, z);
    }
    temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x, z);
    while(_curr_pillar_x - 1 >= 0 && temp_pillar.x * properties.x_mono_dir > x * properties.x_mono_dir) {
        _curr_pillar_x--;
        temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x, z);
    }
    if(_curr_pillar_y + 1 < ny)
        temp_pillar = GetPillarAtZval(_curr_pillar_y + 1, _curr_pillar_x, z);
    while(_curr_pillar_y + 1 < ny && temp_pillar.y * properties.y_mono_dir <= y * properties.y_mono_dir) {
        _curr_pillar_y++;
        temp_pillar = GetPillarAtZval(_curr_pillar_y + 1, _curr_pillar_x, z);
    }
    temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x, z);
    while(_curr_pillar_y - 1 >= 0 && temp_pillar.y * properties.y_mono_dir > y * properties.y_mono_dir) {
        _curr_pillar_y--;
        temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x, z);
    }

    px = _curr_pillar_x;
    py = _curr_pillar_y;

    // 插值计算新的 pillar
    temp_pillar = GetPillarAtZval(_curr_pillar_y, _curr_pillar_x, z);
    temp_pillar1 = GetPillarAtZval(_curr_pillar_y + 1, _curr_pillar_x + 1, z);
    Line pillar_line = GetPillar(_curr_pillar_y, _curr_pillar_x);
    Line pillar_line1 = GetPillar(_curr_pillar_y + 1, _curr_pillar_x + 1); 
    new_pillar.a.x = pillar_line.a.x + (x - temp_pillar.x) / (temp_pillar1.x - temp_pillar.x) * (pillar_line1.a.x - pillar_line.a.x);
    new_pillar.a.y = pillar_line.a.y + (y - temp_pillar.y) / (temp_pillar1.y - temp_pillar.y) * (pillar_line1.a.y - pillar_line.a.y);
    new_pillar.a.z = pillar_line.a.z;
    new_pillar.b.x = pillar_line.b.x + (x - temp_pillar.x) / (temp_pillar1.x - temp_pillar.x) * (pillar_line1.b.x - pillar_line.b.x);
    new_pillar.b.y = pillar_line.b.y + (y - temp_pillar.y) / (temp_pillar1.y - temp_pillar.y) * (pillar_line1.b.y - pillar_line.b.y);
    new_pillar.b.z = pillar_line.b.z;
}

double GRDECLData::InterpolatePillarZvalAtLayer(int layer, Line pillar, INDEX px, INDEX py) {
    return GetCorner(layer, py, px).zval[0][0][0];
}

void GRDECLData::GetNaiveBottomSurface(std::vector<Point3d> &pts, std::vector<INDEX> &bdrs, std::vector<INDEX> &tris) {
    // bottom surface is dense
    pts.reserve((ny + 1) * (nx + 1));
    for (INDEX j = 0; j < ny; j++) {
        for (INDEX i = 0; i < nx; i++) {
            pts.push_back(GetCorner(0, j, i).coord[0][0][0]);
        }
        pts.push_back(GetCorner(0, j, nx - 1).coord[0][0][1]);
    }
    for (INDEX i = 0; i < nx; i++) {
        pts.push_back(GetCorner(0, ny - 1, i).coord[0][1][0]);
    }
    pts.push_back(GetCorner(0, ny - 1, nx - 1).coord[0][1][1]);
    printf("points size: %d\n", pts.size());

    // add boundary lines, including four sides and fault lines

    // add boundary lines
    for (INDEX i = 0; i < nx; i++) {
        bdrs.push_back(i);
        bdrs.push_back(i + 1);
    }

    for (INDEX j = 0; j < ny; j++) {
        bdrs.push_back(j * (nx + 1) + nx);
        bdrs.push_back((j + 1) * (nx + 1) + nx);
    }

    for (INDEX i = 0; i < nx; i++) {
        bdrs.push_back((ny + 1) * (nx + 1) - i - 1);
        bdrs.push_back((ny + 1) * (nx + 1) - i - 2);
    }

    for (INDEX j = 0; j < ny; j++) {
        bdrs.push_back((ny - j) * (nx + 1));
        bdrs.push_back((ny - j - 1) * (nx + 1));
    }

    printf("current boundary line count: %d\n", bdrs.size() / 2);
    // for (int i = 0; i < bdrs.size(); i += 2) {
        // printf("boundary line at (%d %d)\n", bdrs[i], bdrs[i + 1]);
    // }

    // get fault points
    int fault_pts = 0;
    // skip boundary points, only check z = 1 points
    // for (int k = 0; k < nz; k++) { // only consider the bottom surface
        for (int j = 0; j < ny + 1; j++) {
            for (int i = 0; i < nx + 1; i++) {
                if (j > 0 && j < ny && i > 0 && i < nx) { // common inner points
                    if (!isEqual(GetCorner(0, j - 1, i - 1).zval[1][1][1], GetCorner(0, j, i - 1).zval[1][0][1]) ||
                        !isEqual(GetCorner(0, j - 1, i - 1).zval[1][1][1], GetCorner(0, j - 1, i).zval[1][1][0]) ||
                        !isEqual(GetCorner(0, j, i).zval[1][0][0], GetCorner(0, j - 1, i).zval[1][1][0]) ||
                        !isEqual(GetCorner(0, j, i).zval[1][0][0], GetCorner(0, j, i - 1).zval[1][0][1])) {
                        // printf("fault point at (%d, %d)\n", j, i);
                        isFault[j * (nx + 1) + i] = true;
                        fault_pts++;
                        continue;
                    }
                }
                else {
                    // only have inner fault, skip
                }
            }
        }
    // }
    printf("fault points count: %d\n", fault_pts);
    // Warning: the st/ed fault points are not correct

    // add fault lines
    int fault_lines = 0;
    for (int j = 0; j < ny + 1; j++) {
        for (int i = 0; i < nx + 1; i++) {
            if (isFault[j * (nx + 1) + i]) {
                if (isFault[j * (nx + 1) + i + 1]) {
                    // printf("fault line at (%d, %d) and (%d, %d)\n", j, i, j, i + 1);
                    bdrs.push_back(j * (nx + 1) + i);
                    bdrs.push_back(j * (nx + 1) + i + 1);
                    fault_lines++;
                }
                if (isFault[(j + 1) * (nx + 1) + i]) {
                    // printf("fault line at (%d, %d) and (%d, %d)\n", j, i, j + 1, i);
                    bdrs.push_back(j * (nx + 1) + i);
                    bdrs.push_back((j + 1) * (nx + 1) + i);
                    fault_lines++;
                }
            }
        }
    }
    printf("fault line count: %d\n", fault_lines);

    printf("boundary line count: %d\n", bdrs.size() / 2);

    // add tri cells
    tris.reserve(ny * nx * 2 * 3);
    for (INDEX j = 0; j < ny; j++) {
        for (INDEX i = 0; i < nx; i++) {
            // for cell (j, i), add two tri cells, connect (1, 0) and (0, 1) by default
            if (isFault[(j + 1) * (nx + 1) + i] && isFault[j * (nx + 1) + i + 1]) { // connect (0, 0) and (1, 1)
                tris.push_back(j * (nx + 1) + i);
                tris.push_back((j + 1) * (nx + 1) + i + 1);
                tris.push_back((j + 1) * (nx + 1) + i);
                
                tris.push_back(j * (nx + 1) + i);
                tris.push_back((j + 1) * (nx + 1) + i + 1);
                tris.push_back(j * (nx + 1) + i + 1);
            }
            else { // common case? connect (1, 0) and (0, 1)
                tris.push_back(j * (nx + 1) + i);
                tris.push_back(j * (nx + 1) + i + 1);
                tris.push_back((j + 1) * (nx + 1) + i);

                tris.push_back((j + 1) * (nx + 1) + i);
                tris.push_back(j * (nx + 1) + i + 1);
                tris.push_back((j + 1) * (nx + 1) + i + 1);
            }
        }
    }
    printf("tris count: %d\n", tris.size() / 3);

    // Write the VTK file
    std::string bottom_surface_filename = output_path + "/bottom_surface.vtk";
    std::ofstream vtkFile(bottom_surface_filename);
    if (!vtkFile.is_open()) {
        // If the file cannot be opened, print an error message and return
        std::cerr << "Unable to open file: " << bottom_surface_filename << std::endl;
        return;
    }

    // Write the VTK file header
    vtkFile << "# vtk DataFile Version 2.0\n";
    vtkFile << "Grid Legacy VTK File\n";
    vtkFile << "ASCII\n";
    vtkFile << "DATASET UNSTRUCTURED_GRID\n";
    vtkFile << "\n";

    // Write the grid points to the VTK file
    vtkFile << "POINTS " << pts.size() << " float\n";
    for (auto& p : pts) {
        vtkFile << p.x << " " << p.y << " " << p.z << "\n";
    }
    vtkFile << "\n";

    // Write the grid cells to the VTK file
    vtkFile << "CELLS " << tris.size() / 3 << " " << tris.size() + tris.size() / 3 << "\n";
    for (int i = 0; i < tris.size(); i += 3) {
        vtkFile << "3 ";
        vtkFile << tris[i] << " ";
        vtkFile << tris[i + 1] << " ";
        vtkFile << tris[i + 2] << "\n";
    }
    vtkFile << "\n";

    // Write the cell types to the VTK file
    vtkFile << "CELL_TYPES " << tris.size() / 3 << "\n";
    for (int i = 0; i < tris.size() / 3; i++) {
        vtkFile << "5\n";
    }
    vtkFile << "\n";

    // Write point data to the VTK file
    vtkFile << "POINT_DATA " << pts.size() << "\n";
    vtkFile << "SCALARS isFault int 1\n";
    vtkFile << "LOOKUP_TABLE default\n";
    for (int i = 0; i < pts.size(); i++) {
        vtkFile << isFault[i] << "\n";
    }
    vtkFile << "\n";

    // Close the VTK file
    vtkFile.close();

}

void GRDECLData::GetFaultsFromFile(std::vector<Point3d> &fault_points, std::vector<INDEX> &fault_lines, std::vector<INDEX> &fault_quads) {
    double zmin = range.rmin.z;
    double zmax = range.rmax.z;
    std::vector<Point3d> bottom_surface;
    std::vector<Point3d> top_surface;
    bottom_surface.reserve((ny + 1) * (nx + 1));
    top_surface.reserve((ny + 1) * (nx + 1));

    for (INDEX j = 0; j < ny; j++) {
        for (INDEX i = 0; i < nx; i++) {
            bottom_surface.push_back(GetCorner(0, j, i).coord[0][0][0]);
        }
        bottom_surface.push_back(GetCorner(0, j, nx - 1).coord[0][0][1]);
    }
    for (INDEX i = 0; i < nx; i++) {
        bottom_surface.push_back(GetCorner(0, ny - 1, i).coord[0][1][0]);
    }
    bottom_surface.push_back(GetCorner(0, ny - 1, nx - 1).coord[0][1][1]);

    for (INDEX j = 0; j < ny; j++) {
        for (INDEX i = 0; i < nx; i++) {
            top_surface.push_back(GetCorner(nz - 1, j, i).coord[1][0][0]);
        }
        top_surface.push_back(GetCorner(nz - 1, j, nx - 1).coord[1][0][1]);
    }
    for (INDEX i = 0; i < nx; i++) {
        top_surface.push_back(GetCorner(nz - 1, ny - 1, i).coord[1][1][0]);
    }
    top_surface.push_back(GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1]);

    std::vector<std::vector<INDEX>> fault_line_index;
    fault_line_index.resize(faults.size());

    for (int i = 0; i < faults.size(); i++) {
        for (auto line: faults[i].lines) {
            // z value increase from bottom to top by default, and the fault is embedded in the model
            INDEX fault_line = -1; // index for real mesh line
            for (int i = 0; i < line.size() - 1; i++) {
                Point3d &p1 = line[i];
                Point3d &p2 = line[i + 1];
                double tmin = (zmin - p1.z) / (p2.z - p1.z);
                Point3d pmin = {p1.x + tmin * (p2.x - p1.x), p1.y + tmin * (p2.y - p1.y), zmin};
                double tmax = (zmax - p1.z) / (p2.z - p1.z);
                Point3d pmax = {p1.x + tmax * (p2.x - p1.x), p1.y + tmax * (p2.y - p1.y), zmax};
                // printf("fault %s line %d: (%f, %f, %f) to (%f, %f, %f)\n", fault.name.c_str(), i, pmin.x, pmin.y, pmin.z, pmax.x, pmax.y, pmax.z);

                INDEX line_index;
                double line_dist;
                findClosestLine(bottom_surface, top_surface, pmin, pmax, line_index, line_dist);
                printf("line index: (%d, %d), dist: %f\n", line_index / (nx + 1), line_index % (nx + 1), line_dist);

                if (fault_line == -1 || fault_line == line_index) {
                    fault_line = line_index;
                }
                else { // maybe error, but never happened in case gridA/B
                    printf("Error: single fault line has multiple line results\n");
                    break;
                }
            }
            fault_line_index[i].push_back(fault_line);
            isFault[fault_line] = true; // fault surface point
        }
    }

    // std::vector<std::vector<INDEX>> fault_surfaces(fault_line_index); // real output

    // calculate fault intersection points
    if (fault_line_index.size() > 1) {
        for (int i = 0; i < fault_line_index.size(); i++) {
            for (int j = i + 1; j < fault_line_index.size(); j++) {
                for (int m = 0; m < fault_line_index[i].size() - 1; m++) {
                    for (int n = 0; n < fault_line_index[j].size() - 1; n++) {
                        // TODO: using function instead
                        Point2d pmin; // bottom surface intersection point
                        Point2d pmin1 = {bottom_surface[fault_line_index[i][m]].x, bottom_surface[fault_line_index[i][m]].y};
                        Point2d pmin2 = {bottom_surface[fault_line_index[i][m + 1]].x, bottom_surface[fault_line_index[i][m + 1]].y};
                        Point2d pmin3 = {bottom_surface[fault_line_index[j][n]].x, bottom_surface[fault_line_index[j][n]].y};
                        Point2d pmin4 = {bottom_surface[fault_line_index[j][n + 1]].x, bottom_surface[fault_line_index[j][n + 1]].y};

                        Point2d pmax; // top surface intersection point
                        Point2d pmax1 = {top_surface[fault_line_index[i][m]].x, top_surface[fault_line_index[i][m]].y};
                        Point2d pmax2 = {top_surface[fault_line_index[i][m + 1]].x, top_surface[fault_line_index[i][m + 1]].y};
                        Point2d pmax3 = {top_surface[fault_line_index[j][n]].x, top_surface[fault_line_index[j][n]].y};
                        Point2d pmax4 = {top_surface[fault_line_index[j][n + 1]].x, top_surface[fault_line_index[j][n + 1]].y};

                        if (segmentsIntersect(pmin1, pmin2, pmin3, pmin4) && segmentsIntersect(pmax1, pmax2, pmax3, pmax4)) {
                            if (lineIntersection(pmin1, pmin2, pmin3, pmin4, pmin)) {
                                printf("bottom surface intersection point: (%f, %f)\n", pmin.x, pmin.y);
                            }

                            if (lineIntersection(pmax1, pmax2, pmax3, pmax4, pmax)) {
                                printf("top surface intersection point: (%f, %f)\n", pmax.x, pmax.y);
                            }
                        }
                        else {
                            continue;
                        }

                        INDEX line_index;
                        double line_dist;
                        findClosestLine(bottom_surface, top_surface, {pmin.x, pmin.y, zmin}, {pmax.x, pmax.y, zmax}, line_index, line_dist);
                        printf("line index: (%d, %d), dist: %f\n", line_index / (nx + 1), line_index % (nx + 1), line_dist);
                        // Warning: vector length may change, but work in grid b
                        fault_line_index[i].insert(fault_line_index[i].begin() + m + 1, line_index); m++;
                        fault_line_index[j].insert(fault_line_index[j].begin() + n + 1, line_index); n++;
                    }
                }
            }
        }
    }

    // TODO: support grid b
    // add fault point
    const INDEX layer_z = nz + 1;
    for (int i = 0; i < fault_line_index.size(); i++) {
        for (int j = 0; j < fault_line_index[i].size(); j++) {
            // 0 and size() - 1 are boundary dense points
            INDEX line_y = fault_line_index[i][j] / (nx + 1);
            INDEX line_x = fault_line_index[i][j] % (nx + 1);
            printf("fault %d line %d: (%d, %d)\n", i, j, line_y, line_x);

            if (j == 0 || j == fault_line_index[i].size() - 1) {
                // inner the surface
                for (int k = 0; k < layer_z; k++) {
                    Point3d p = GetCorner(k, line_y, line_x).coord[0][0][0];
                    fault_points.push_back(p);
                }
            }
            else { // not dense point

            }
        }
    }
}

void GRDECLData::GetCornerOfPoint(INDEX &z, INDEX &y, INDEX &x,
                                  const Point3d &p, INDEX prev_y, INDEX prev_x, std::mt19937 &gen) {
    double zval = p.z;
    x = prev_x;
    y = prev_y;

    while (true) {
        // 取当前单元4个角点（按顺序）
        // example: 2 - - - 3
        //          |       |
        //          |       |
        //          0 - - - 1
        Point2d p1 = GetPointFromPillar(y,     x,     zval);
        Point2d p2 = GetPointFromPillar(y,     x + 1, zval);
        Point2d p3 = GetPointFromPillar(y + 1, x,     zval);
        Point2d p4 = GetPointFromPillar(y + 1, x + 1, zval);

        auto convex = checkQuadrilateralConvexity(p1, p2, p3, p4);

        if (convex.first) { // 四边形凸（正常情况）
            // 计算四角中心作为参考点
            Point2d p_ref = (p1 + p2 + p3 + p4) / 4;
            double ref_cross1 = crossProduct(p1, p2, p_ref);
            double ref_cross2 = crossProduct(p2, p4, p_ref);
            double ref_cross3 = crossProduct(p4, p3, p_ref);
            double ref_cross4 = crossProduct(p3, p1, p_ref);

            // 检查参考点是否在单元内
            if (!(ref_cross1 < 0 && ref_cross2 < 0 && ref_cross3 < 0 && ref_cross4 < 0) &&
                !(ref_cross1 > 0 && ref_cross2 > 0 && ref_cross3 > 0 && ref_cross4 > 0)) {
                printf("Error: ref point is not in the grid\n");
            }

            // 计算 p 与各边的叉积
            double cross1 = crossProduct(p1, p2, {p.x, p.y});
            double cross2 = crossProduct(p2, p4, {p.x, p.y});
            double cross3 = crossProduct(p4, p3, {p.x, p.y});
            double cross4 = crossProduct(p3, p1, {p.x, p.y});

            if ((cross1 < 0 && cross2 < 0 && cross3 < 0 && cross4 < 0) ||
                (cross1 > 0 && cross2 > 0 && cross3 > 0 && cross4 > 0)) {
                // p 位于当前四边形内，接下来确定 z 方向索引
                double z_bottom;
                double z_top;

                for (int i = 0; i < nz; i++) {
                    z_bottom = (GetCorner(i, y, x).zval[0][0][0] + GetCorner(i, y, x).zval[0][0][1] +
                                GetCorner(i, y, x).zval[0][1][0] + GetCorner(i, y, x).zval[0][1][1]) / 4;
                    z_top = (GetCorner(i, y, x).zval[1][0][0] + GetCorner(i, y, x).zval[1][0][1] +
                                GetCorner(i, y, x).zval[1][1][0] + GetCorner(i, y, x).zval[1][1][1]) / 4;

                    if ((zval - z_bottom) * (zval - z_top) <= 0) {
                        z = i;
                        break;
                    }
                }
                // if (z == -1) {
                    // printf("zval: %f, unreachable branch: (%d, %d)\n", zval, y, x);
                // }
                break; // 找到所属单元，退出循环
            }
            // printf("convex branch: (%d, %d)\n", y, x);

            // 如果 p 不在当前单元内，则根据 p 与各边的相对位置调整单元索引
            if (cross1 * ref_cross1 < 0 && cross3 * ref_cross3 > 0 && y > 0) {
                y--;
            }
            else if (cross1 * ref_cross1 > 0 && cross3 * ref_cross3 < 0 && y < ny - 1) {
                y++;
            }

            if (cross2 * ref_cross2 < 0 && cross4 * ref_cross4 > 0 && x < nx - 1) {
                x++;
            }
            else if (cross2 * ref_cross2 > 0 && cross4 * ref_cross4 < 0 && x > 0) {
                x--;
            }
        }
        else { 
            // 四边形为凹：我们将其沿凹角对角线切分为两个三角形
            // 注意：convex.second 返回 1~4，对应 p1~p4中“中间点”，我们用其减 1 得到凹角索引
            // 注意：注意对角线顶点编号
            Point2d arr[4] = { p1, p2, p3, p4 };
            int concaveIndex = convex.second - 1; // 取出凹角顶点索引（0~3）
            int oppIndex   = 3 - concaveIndex; // 对应凹角对面的顶点
            int nextIndex  = concaveIndex ^ 1; // 对应凹角的下一个顶点
            int prevIndex  = concaveIndex ^ 2; // 对应凹角的前一个顶点

            // 分别构成两个三角形：
            // 三角形1： (arr[concaveIndex], arr[nextIndex], arr[oppIndex])
            // 三角形2： (arr[concaveIndex], arr[oppIndex], arr[prevIndex])

            double tri11 = crossProduct(arr[concaveIndex], arr[nextIndex], {p.x, p.y});
            double tri12 = crossProduct(arr[nextIndex], arr[oppIndex], {p.x, p.y});
            double tri13 = crossProduct(arr[oppIndex], arr[concaveIndex], {p.x, p.y}); // new line

            double tri21 = crossProduct(arr[concaveIndex], arr[oppIndex], {p.x, p.y}); // new line
            double tri22 = crossProduct(arr[oppIndex], arr[prevIndex], {p.x, p.y});
            double tri23 = crossProduct(arr[prevIndex], arr[concaveIndex], {p.x, p.y});

            bool inTri1 = (tri11 < 0 && tri12 < 0 && tri13 < 0) || (tri11 > 0 && tri12 > 0 && tri13 > 0);
            bool inTri2 = (tri21 < 0 && tri22 < 0 && tri23 < 0) || (tri21 > 0 && tri22 > 0 && tri23 > 0);

            if (inTri1 || inTri2) {
                // 如果 p 落在其中一个三角形内，则认为 p 在当前单元内
                // p 位于当前四边形内，接下来确定 z 方向索引
                double z_bottom;
                double z_top;

                for (int i = 0; i < nz; i++) {
                    z_bottom = (GetCorner(i, y, x).zval[0][0][0] + GetCorner(i, y, x).zval[0][0][1] +
                                GetCorner(i, y, x).zval[0][1][0] + GetCorner(i, y, x).zval[0][1][1]) / 4;
                    z_top = (GetCorner(i, y, x).zval[1][0][0] + GetCorner(i, y, x).zval[1][0][1] +
                                GetCorner(i, y, x).zval[1][1][0] + GetCorner(i, y, x).zval[1][1][1]) / 4;

                    if ((zval - z_bottom) * (zval - z_top) <= 0) {
                        z = i;
                        break;
                    }
                }
                // if (z == -1) {
                    // printf("zval: %f, unreachable branch: (%d, %d)\n", zval, y, x);
                // }
                break; // 找到所属单元，退出循环
            }
            else {
                // printf("concave branch: (%d, %d)\n", y, x);
                // hard to determine the next point, use random walk
                x = x + randomMinusOneZeroOne(gen);
                y = y + randomMinusOneZeroOne(gen);

            }
        } // end concave branch
    } // end while

    // printf("final corner index: (%d, %d, %d)\n", z, y, x);
}
